In [222]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
#sns.set()

In [223]:
import numpy as np
import h5py
from chainconsumer import ChainConsumer
from glob import glob
from os import path
from copy import deepcopy
from collections import OrderedDict

In [224]:
! ls -ltr ~/des/PearceMCMC/UniverseMachine*HOD.hdf5


-rw-r--r-- 1 swmclau2 des 324473989 Feb 26 01:48 /u/ki/swmclau2/des/PearceMCMC/UniverseMachine_wp_ds_rmin_1.0_HOD.hdf5
-rw-r--r-- 1 swmclau2 des 373731747 Feb 26 01:48 /u/ki/swmclau2/des/PearceMCMC/UniverseMachine_wp_ds_rmin_2.0_HOD.hdf5
-rw-r--r-- 1 swmclau2 des 321548554 Feb 26 01:48 /u/ki/swmclau2/des/PearceMCMC/UniverseMachine_wp_ds_rmin_0.5_HOD.hdf5
-rw-r--r-- 1 swmclau2 des     19664 Feb 26 01:48 /u/ki/swmclau2/des/PearceMCMC/UniverseMachine_wp_ds_rmin_5.0_HOD.hdf5
-rw-r--r-- 1 swmclau2 des   1613847 Feb 26 09:45 /u/ki/swmclau2/des/PearceMCMC/UniverseMachine_wp_ds_rmin_None_HOD.hdf5

In [225]:
fnames = sorted(glob('/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_*_CAB.hdf5'))

Assume all chains run with same true values, only thing that differs is the rmin.


In [226]:
fnames


Out[226]:
['/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_0.5_CAB.hdf5',
 '/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_1.0_CAB.hdf5',
 '/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_2.0_CAB.hdf5',
 '/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_None_CAB.hdf5']

In [227]:
fnames = [ i for i in reversed(fnames[:-1])]
_fnames = [] _fnames.extend(reversed(fnames[:-1]) ) _fnames.append(fnames[-1]) fnames = _fnames

In [228]:
fnames


Out[228]:
['/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_2.0_CAB.hdf5',
 '/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_1.0_CAB.hdf5',
 '/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_0.5_CAB.hdf5']

In [229]:
model = path.basename(fnames[0]).split('_')[-1].split('.')[0]
#model = path.basename(fnames[0]).split('_')[-2]#.split('.')[0]

In [230]:
model


Out[230]:
'CAB'

In [231]:
f = h5py.File(fnames[-1], 'r')

In [232]:
f.attrs.keys()


Out[232]:
[u'chain_fixed_params',
 u'cov',
 u'dlogz',
 u'emu_cov_fname',
 u'emu_hps',
 u'emu_type',
 u'fixed_params',
 u'mcmc_type',
 u'nburn',
 u'nlive',
 u'nsteps',
 u'nwalkers',
 u'obs',
 u'param_names',
 u'seed',
 u'sim',
 u'training_file',
 u'true_cov_fname',
 u'true_data_fname']

In [233]:
n_walkers = f.attrs['nwalkers']

In [234]:
chain_pnames = f.attrs['param_names']

In [235]:
#sim_info = eval(f.attrs['sim'])
sim_info = f.attrs['sim']

In [236]:
#gal_type = eval(f.attrs['sim'])['gal_type']
gal_type = 'SHAM'

In [237]:
f.close()

In [238]:
chain_pnames


Out[238]:
array(['ombh2', 'omch2', 'w0', 'ns', 'H0', 'Neff', 'sigma8',
       'mean_occupation_centrals_assembias_param1',
       'mean_occupation_satellites_assembias_slope1', 'logM0',
       'sigma_logM', 'mean_occupation_satellites_assembias_param1',
       'mean_occupation_centrals_assembias_slope1', 'logM1', 'alpha',
       'conc_gal_bias'], dtype='|S43')

In [239]:
param_name_dict = {'ombh2': r'$\Omega_b h^2$', 'omch2': r'$\Omega_c h^2$','w0': r'$w_0$','ns': r'$n_s$', \
                   'sigma8': r'$\sigma_8$', 'ln10As':r'$\ln 10^{10} A_s$', 'H0': r'$H_0$','Neff': r'$N_{eff}$',\
               'mean_occupation_centrals_assembias_corr1': r'$\rho_{cen}$',\
                 'mean_occupation_satellites_assembias_corr1':r'$\rho_{sat}$',\
                   'mean_occupation_centrals_assembias_param1': r'$\mathcal{A}_{cen}$',\
                 'mean_occupation_satellites_assembias_param1':r'$\mathcal{A}_{sat}$',\
                   'mean_occupation_centrals_assembias_slope1': r'$\mathcal{B}_{cen}$',\
                 'mean_occupation_satellites_assembias_slope1':r'$\mathcal{B}_{sat}$',\
                   'logM1': r'$\log(M_1)$','logM0': r'$\log(M_0)$','sigma_logM': r'$\sigma_{\log M }$',
                   'conc_gal_bias': r'$\eta$', 'alpha':r'$\alpha$' }

In [240]:
bounds_dict = {'H0': (61.69472, 74.76751999999999),
 'Neff': (2.62125, 4.27875),
 'alpha': (0.7, 1.3),
 'conc_gal_bias': (0.5, 2.0),
 'sigma8': (0.65, 1.0), # TODO update
 'logM0': (12.6, 13.6),
 'logM1': (13.7, 14.7),
 'ns': (0.9278462, 0.9974495999999999),
 'ombh2': (0.02066455, 0.02371239),
 'omch2': (0.1012181, 0.13177679999999997),
 'sigma_logM': (0.05, 0.5),
 'ln10As': (3.0, 3.1),
 'w0': (-1.399921, -0.5658486),
 'mean_occupation_centrals_assembias_corr1': (-1.0, 1.0),
 'mean_occupation_satellites_assembias_corr1': (-1.0, 1.0),
 'mean_occupation_centrals_assembias_param1': (-1.0, 1.0),
 'mean_occupation_satellites_assembias_param1': (-1.0, 1.0),
 'mean_occupation_centrals_assembias_slope1': (-3.0, 3.0),
 'mean_occupation_satellites_assembias_slope1': (-3.0, 3.0)}

In [241]:
hod_param_names = []
cosmo_param_names = []

hod_bounds = []
cosmo_bounds = []
cosmo_names = set(['ombh2', 'omch2', 'w0', 'ns', 'sigma8', 'H0', 'Neff', 'ln10As'])
for pname in chain_pnames:
    if pname in cosmo_names:
        cosmo_param_names.append(param_name_dict[pname])
        cosmo_bounds.append(bounds_dict[pname])
    else:
        hod_param_names.append(param_name_dict[pname])
        hod_bounds.append(bounds_dict[pname])
        
param_names = deepcopy(cosmo_param_names)
param_names.extend(hod_param_names)

bounds = deepcopy(cosmo_bounds)
bounds.extend(hod_bounds)

In [242]:
c = ChainConsumer()

In [243]:
n_burn = 1000

for fname in fnames:
#for fname in [fnames[0], fnames[3]]:
    try:
        f = h5py.File(fname, 'r')
        chain = f['chain'][n_burn*n_walkers:]
        chain = chain[np.all(chain!=0.0, axis = 1), :]
    except IOError:
        print 'Error loading', fname
        continue
    
    #print chain.shape
    chain = chain.reshape((-1, n_walkers, chain.shape[1]))
    #print chain.shape
    chain = chain.reshape((-1, chain.shape[2]), order = 'F')
    #chain = chain[:, 0, :]
    if chain.shape[0] == 0:
        print 'Insufficient samples in', fname
        continue
    
    #name = path.basename(fname).split('_')[-1].split('.')[0]
    name = path.basename(fname).split('_')[4]
    #name = path.basename(fname).split('_')[3]

    if name == "None":
        name = "0.1"
        
    name = r"$r_{min} = %s$"%name
    
    hod_param_names = []
    cosmo_param_names = []
    chain_pnames = f.attrs['param_names']
    cosmo_names = set(['ombh2', 'omch2', 'w0', 'ns', 'sigma8', 'H0', 'Neff'])
    for pname in chain_pnames:
        if pname in cosmo_names:
            cosmo_param_names.append(param_name_dict[pname])
        else:
            hod_param_names.append(param_name_dict[pname])

    param_names = deepcopy(cosmo_param_names)
    param_names.extend(hod_param_names)
    print fname
    print chain.shape, chain.shape[0]/n_walkers, len(param_names)
    c.add_chain(chain, parameters=param_names, name = name, walkers = n_walkers)
    f.close()


Insufficient samples in /u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_2.0_CAB.hdf5
/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_1.0_CAB.hdf5
(1163000, 16) 4652 16
/u/ki/swmclau2/des/PearceMCMC/ShuffledSHAM_wp_ds_rmin_0.5_CAB.hdf5
(2395250, 16) 9581 16
for start in np.arange(0.0, 3.0, 0.2): palette = sns.cubehelix_palette(len(fnames), start = start, rot = 0.2, gamma = 0.7) print start sns.palplot(palette) plt.show()

In [244]:
if model == 'HOD':
    #palette = sns.cubehelix_palette(len(fnames), start = 2.6, rot = 0.1, gamma = 0.7)
    palette = sns.cubehelix_palette(len(fnames), start = 2.6, rot = 0.1, gamma = 1.2, reverse=False)

elif model == 'HSAB':
    palette = sns.cubehelix_palette(len(fnames), start = 2.0, rot = 0.1, gamma = 0.7)
    #palette = sns.cubehelix_palette(len(fnames), start = 1.2, rot = 0.6, gamma = 1.7)
elif model == 'CAB':
    palette = sns.cubehelix_palette(len(fnames), start = 0.6, rot = 0.1, gamma = 1.2)
elif model == 'CorrAB':
    #palette = sns.cubehelix_palette(len(fnames), start = 1.2, rot = 0.1, gamma = 1.2)
    palette = sns.cubehelix_palette(len(fnames), start = 1.2, rot = 0.1, gamma = 2.0)

sns.palplot(palette)



In [245]:
def color_to_hex(color):
    return np.array(color)#*255

In [246]:
c.configure(colors = [color_to_hex(p) for p in palette], shade = True, shade_alpha=0.2, shade_gradient=1.0)


Out[246]:
<chainconsumer.chainconsumer.ChainConsumer at 0x7fdaa4fe8d90>

In [247]:
n_params = chain.shape[1] if len(chain.shape) > 1 else 1

In [248]:
hod_idxs = np.array(range(len(cosmo_param_names), len(cosmo_param_names)+len(hod_param_names)-1)) #skip eta
cosmo_idxs = np.array(range(len(cosmo_param_names)))

In [249]:
cosmo_param_names


Out[249]:
['$\\Omega_b h^2$',
 '$\\Omega_c h^2$',
 '$w_0$',
 '$n_s$',
 '$H_0$',
 '$N_{eff}$',
 '$\\sigma_8$']

In [250]:
hod_param_names


Out[250]:
['$\\mathcal{A}_{cen}$',
 '$\\mathcal{B}_{sat}$',
 '$\\log(M_0)$',
 '$\\sigma_{\\log M }$',
 '$\\mathcal{A}_{sat}$',
 '$\\mathcal{B}_{cen}$',
 '$\\log(M_1)$',
 '$\\alpha$',
 '$\\eta$']

In [251]:
if chain.shape[1] == 7:
    cosmo_chain = chain
    
elif chain.shape[1] == 5:
    hod_chain = chain
else:
    hod_chain = chain[:,hod_idxs]
    #cosmo_chain = chain[:,cosmo_idxs]

In [252]:
gal_type = 'SHAM'

In [253]:
if gal_type == 'HOD':
    from pearce.mocks import cat_dict
    #cosmo_params = {'simname': sim_info['simname'], 'boxno': sim_info['sim_hps']['boxno'],\
    #                'realization': sim_info['sim_hps']['realization'], 'scale_factors':[sim_info['scale_factor']],\
    #                'system': 'ki-ls'}
    cosmo_params = {'simname': 'testbox', 'boxno': 1,\
                    'realization':0, 'scale_factors':[1.0],\
                    'system': 'ki-ls'}
    cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!
    
    cpv = cat._get_cosmo_param_names_vals()

    cat_val_dict =  {key: val for key, val in zip(cpv[0], cpv[1])}
    cosmo_true_vals = [cat_val_dict[pn] for pn in chain_pnames if pn in cat_val_dict]

    hod_params = {'alpha': 1.083, 'conc_gal_bias': 1.0, 'logM0': 13.2,'logM1': 14.2, 'sigma_logM': 0.2}#sim_info['hod_params']
    hod_true_vals = [hod_params[key] for key in chain_pnames if key in hod_params]
    
    if model!= 'HOD':
        hod_true_vals.extend([0, 0])
        if model == 'CAB':
            hod_true_vals.extend([0,0])
else: #SHAM, etc on MDPL2
 
    # multidark planck2
    #cosmo_true_vals = np.array([0.0223, 0.1188, -1, 0.9667, 3.047, \
    #                            0.6774*100, 3.046]) #darksky
    
    h = 0.6777
    #cosmo_pnames = ['ombh2', 'omch2', 'w0', 'ns', 'H0', 'Neff', 'sigma8']
    cosmo_true_vals = np.array([0.048206*h**2, 0.307115*h**2- 0.048206*h**2,\
                                -1, 0.9667, \
                                h*100, 3.046, 0.8228]) #mdpl2

    hod_true_vals = np.array([np.inf for p in hod_param_names])
c.plotter.plot_walks(parameters = [cosmo_param_names[1]],\ truth=[cosmo_true_vals[1]]);#, convolve = 10);

In [254]:
hod_param_names


Out[254]:
['$\\mathcal{A}_{cen}$',
 '$\\mathcal{B}_{sat}$',
 '$\\log(M_0)$',
 '$\\sigma_{\\log M }$',
 '$\\mathcal{A}_{sat}$',
 '$\\mathcal{B}_{cen}$',
 '$\\log(M_1)$',
 '$\\alpha$',
 '$\\eta$']
# Enforce a consistent order for the plots #if gal_type == 'HOD': plot_hod_param_names = ['$\\sigma_{\\log M }$', '$\\log(M_0)$', '$\\log(M_1)$', '$\\alpha$', '$\\eta$'] plot_hod_true_vals = [0.2, 13.2, 14.2, 1.083, 1.0] plot_hod_bounds = [(0.05, 0.5), (12.6, 13.6), (13.7, 14.7), (0.7, 1.3), (0.5, 2.0)] ab_param_names = [r'$\mathcal{A}_{%s}$', r'$\mathcal{B}_{%s}$', r'$\rho_{%s}$'] ab_true_vals = [0.0, np.inf,0.0] for abpn, hodtv in zip(ab_param_names, ab_true_vals): if abpn%'cen' in hod_param_names: # this pname is in the model for gal_type in ['cen', 'sat']: plot_hod_param_names.append(abpn%gal_type) if gal_type == 'cen': plot_hod_true_vals.append(1.0) else: plot_hod_true_vals.append(-1.0) plot_hod_bounds.append((-1.0, 1.0)) #for phpn in plot_hod_param_names: # i = hod_param_names.index(phpn) # print phpn, i, hod_true_vals[i] # plot_hod_true_vals.append(hod_true_vals[i]) # plot_hod_bounds.append(hod_bounds[i])

In [255]:
# Enforce a consistent order for the plots
#if gal_type == 'HOD':
    
plot_hod_param_names = ['$\\sigma_{\\log M }$', '$\\log(M_0)$', '$\\log(M_1)$', '$\\alpha$', '$\\eta$']
plot_hod_true_vals = [0.2, 13.2, 14.2, 1.083, 1.0]
#plot_hod_true_vals = [0.44415644, 11.41003864, 14.56088772, 1.03887697, 1.0]

plot_hod_bounds = [(0.05, 0.5), (12.6, 13.6), (13.7, 14.7), (0.7, 1.3), (0.5, 2.0)]

ab_param_names = [r'$\mathcal{A}_{%s}$', r'$\mathcal{B}_{%s}$', r'$\rho_{%s}$']
ab_true_vals = [0.0, np.inf,0.0]
for abpn, hodtv in zip(ab_param_names, ab_true_vals):
    if abpn%'cen' in hod_param_names: # this pname is in the model
        if 'B' in abpn:
            plot_hod_true_vals.extend([np.inf, np.inf])
            plot_hod_param_names.extend([abpn%'cen', abpn%'sat'])

            continue
        for gal_type in ['cen', 'sat']:
            plot_hod_param_names.append(abpn%gal_type)#

            if gal_type == 'cen':
                plot_hod_true_vals.append(1.0)
            else:
                plot_hod_true_vals.append(-1.0)

            plot_hod_bounds.append((-1.0, 1.0))

#for phpn in plot_hod_param_names:
#    i = hod_param_names.index(phpn)
#    print phpn, i, hod_true_vals[i]
#    plot_hod_true_vals.append(hod_true_vals[i])
#    plot_hod_bounds.append(hod_bounds[i])
plot_hod_true_vals

In [256]:
fig = c.plotter.plot(figsize='GROW', parameters = plot_hod_param_names, truth=plot_hod_true_vals)#, extents=hod_bounds) 
#plt.suptitle('HSAB')
fig.show()


WARNING:chainconsumer.analysis:Parameter $\log(M_0)$ in chain $r_{min} = 1.0$ is not constrained
WARNING:chainconsumer.analysis:Parameter $\alpha$ in chain $r_{min} = 0.5$ is not constrained
WARNING:chainconsumer.analysis:Parameter $\eta$ in chain $r_{min} = 0.5$ is not constrained
WARNING:chainconsumer.analysis:Parameter $\mathcal{B}_{sat}$ in chain $r_{min} = 1.0$ is not constrained
WARNING:chainconsumer.analysis:Parameter $\mathcal{B}_{sat}$ in chain $r_{min} = 0.5$ is not constrained
/u/ki/swmclau2/.conda/envs/hodemulator/lib/python2.7/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
fig = c.plotter.plot(figsize='GROW',parameters = plot_hod_param_names, truth=plot_hod_true_vals, extents=plot_hod_bounds) #plt.suptitle('HSAB') fig.show()
fig = c.plotter.plot_summary(parameters = plot_hod_param_names, truth=plot_hod_true_vals,\ figsize=2, errorbar=True) #fig.suptitle(model) # + galtype fig.show()

In [257]:
cosmo_true_vals


Out[257]:
array([ 2.21399210e-02,  1.18911024e-01, -1.00000000e+00,  9.66700000e-01,
        6.77700000e+01,  3.04600000e+00,  8.22800000e-01])
cosmo_chain.mean(axis = 0)

In [258]:
fig = c.plotter.plot(figsize='PAGE', parameters = cosmo_param_names, truth=cosmo_true_vals, extents=cosmo_bounds) 
fig.show()


WARNING:chainconsumer.analysis:Parameter $\Omega_b h^2$ in chain $r_{min} = 0.5$ is not constrained
WARNING:chainconsumer.analysis:Parameter $n_s$ in chain $r_{min} = 0.5$ is not constrained
WARNING:chainconsumer.analysis:Parameter $N_{eff}$ in chain $r_{min} = 1.0$ is not constrained
WARNING:chainconsumer.analysis:Parameter $N_{eff}$ in chain $r_{min} = 0.5$ is not constrained
/u/ki/swmclau2/.conda/envs/hodemulator/lib/python2.7/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "

In [259]:
from copy import deepcopy
true_vals  = list(deepcopy(cosmo_true_vals))
true_vals.extend(plot_hod_true_vals)

In [260]:
plot_param_names = list(deepcopy(cosmo_param_names))
plot_param_names.extend(plot_hod_param_names)

In [261]:
plot_param_bounds = list(deepcopy(cosmo_bounds))
plot_param_bounds.extend(plot_hod_bounds)

In [262]:
len(plot_param_names)


Out[262]:
16

In [263]:
# TODO can choose these by galaxy model 
# TODO combine multiple chains? 
#summary_idxs = [1, 4]
summary_idxs = [1, 4,6]
if len(plot_param_names) > 12:
    summary_idxs.extend([12,13])
    #summary_idxs.extend([11,12])
    #pass


summary_pnames = [plot_param_names[i] for i in summary_idxs]
summary_truths = [true_vals[i] for i in summary_idxs]
summary_bounds = [plot_param_bounds[i] for i in summary_idxs]

In [264]:
fig = c.plotter.plot(figsize='PAGE', parameters =summary_pnames\
                                    , truth=summary_truths, extents=summary_bounds) 
fig.show()


/u/ki/swmclau2/.conda/envs/hodemulator/lib/python2.7/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "

In [265]:
fig = c.plotter.plot_summary(parameters = summary_pnames, truth=summary_truths, extents=summary_bounds,\
                             figsize=2, errorbar=True) 
#fig.suptitle(model) # + galtype
fig.show()


/u/ki/swmclau2/.conda/envs/hodemulator/lib/python2.7/site-packages/matplotlib/figure.py:403: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "
fig = c.plotter.plot(figsize='GROW', parameters = plot_param_names, truth=true_vals)#, extents=bounds) fig.show()
cosmo_true_vals
cat.cosmology
import pyccl as ccl
cosmo_param_names
ShuffledSHAMdef compute_s8(sample): sample = sample.astype(float) h = sample[5]/100 Omb = sample[0]/h**2 Omc = sample[1]/h**2 w0 = sample[2] n_s = sample[3] A_s = np.exp(sample[4])*(10**(-10)) Neff = sample[6] #print cosmo_param_names #print Omb, Omc, h, A_s, n_s, Neff, w0 cosmo = ccl.Cosmology(Omega_c = Omc, Omega_b = Omb, h=h, A_s = A_s, n_s = n_s, Neff = Neff, w0=w0) sigma_8 = ccl.sigma8(cosmo) Om = (Omb+Omc) s8 = sigma_8*np.sqrt(Om/0.3) return np.array([Om, sigma_8, s8])
s8_chain = [] for i, sample in enumerate(cosmo_chain): print i, s8_sample = compute_s8(cosmo_chain[0,:]) s8_chain.append(s8_sample)

In [ ]:


In [ ]:


In [ ]: